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Abstract 

The Abelian Sandpile Model (ASM) is a paradigm of self-organized criticality (SOC) which 
is related to c = — 2 conformal field theory. The conformal fields corresponding to some height 
clusters have been suggested before. Here we derive the first corrections to such fields, in a 
field theoretical approach, when the lattice parameter is non-vanishing and consider them in 
the presence of a boundary. 
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1 Introduction 

Self-organized criticality is believed to be the underlying reason for the scaling laws seen in a 
number of natural phenomena PQ. Dynamics of a self-organized critical system is such that the 
system naturally approaches its critical state and exhibits long range orders and scaling laws. This 
means that, unlike most of statistical models, without any fine-tuning of some parameters such as 
temperature, the system reaches its critical point. 

The concept was first introduced by Bak, Tang and Wiesenfeld [2j. In their paper they proposed 
the Abelian Sandpile Model (ASM) as a model of self-organized criticality. Since then many different 
models exhibiting the same phenomenon have been developed, but still ASM is the simplest, most 
studied model, in which many analytical results has been derived. For a good review see ref. [3]. 

Many exact results are derived in this theory. The first analytical calculation, which paved 
the road for other analytical results, was done by Dhar [5]. Probabilities of some specific clusters, 
known as Weakly Allowed Clusters (WAC's), was then computed [6]. The simplest of these clusters 
is one-site height one cluster. The probabilities of other one-site clusters with height above 1 was 
computed in ref. [7]. Among other analytical results one can mention the results on boundary 
correlations of height variables and effect of boundary conditions [HJ [TU1 [TT] I12| . on presence of 
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dissipation in the model [10l [HI 13] , on field theoretical approaches [U [TUl 131 IS], on finite size 
corrections [12j [6] and many other results [3]. 

On the other hand, the model has been related to some other lattice models such as spanning 
trees, which correspond to the well-known c = — 2 conformal model. Mahieu and Ruelle [3] related 
ASM to logarithmic conformal field theory through introducing scaling fields corresponding to 
the Weakly Allowed Clusters (WACs), such as one site height-1 cluster. This was done through 
comparison of correlation functions in the two models, and the fields corresponding to a list of 
simplest WACs was obtained. Later, Jeng [9j introduced an elegant method to calculate such 
fields for arbitrary WAC and showed that for all of these clusters the the scaling dimension is 
two. Meanwhile, a more direct way to show this correspondence was developed in ref. [13] • The 
benefit of this new method is that it comes from an action and is not found merely by comparing 
correlation functions. This allows some further investigations that was not possible before. The 
fields associated with other one-site clusters are also derived [131 [TJl [TU] . 

However, most of these results are obtained in the thermodynamic limit; i.e., it is supposed that 
the size of the system, L, is very large and the lattice spacing, o, is ignorable. Though some results 
are obtained in situations where these conditions are not satisfied, we would like to use the method 
introduced in ref. [13] to obtain scaling fields when we are far from thermodynamic limit. To loose 
the first assumption {L — > oo) one should consider the finite size effects and to see the effects of 
loosing the second condition, one may consider how the scaling fields depend on the lattice spacing. 
Because of discrete nature of the model, it is important to investigate this problem. 

In this paper we first review the method introduced in [14] to derive scaling fields and then 
using the method, we derive the higher order corrections to the first order of lattice spacing. In 
the end we discuss the effect of boundary on the field derived using the method of [13] • 

2 Scaling Fields 

As mentioned before, a correspondence between ASM and other well-known statistical models has 
been established, the most significant of which are the connection to q — > limit of the q-state 
Potts model [IS], spanning trees [B] and dense polymers (T8J. All of these models display conformal 
symmetry at their critical points and the proposed CFT corresponding to them is c = —2 model, 
which is a logarithmic CFT [19]. Therefore, it is reasonable to seek a straight way to connect ASM 
to this conformal field theory. One of the first attempts was done by Mahieu and Ruelle. Through 
some arguments of locality and scaling dimension and comparison of height correlations of ASM 
with correlation functions of different fields in c = —2 model, they noted that certain fields can be 
associated to a number of clusters; namely, WACs [4J. Though it was a great step forward, the 
method still had the shortcoming that it was based only on comparison of correlation functions. 
Inspired by their results and a method suggested by Ivashkevich to relate dense polymers and CFT 
|16] - later elaborated for all trees and forests by Caracciolo et aZ. [17j— Moghimi-Araghi, Rajabpour 
and Rouhani solidly confirmed the correspondence between c = —2 and ASM height correlations 
|14j . They deduced the results of Mahieu and Ruelle [4] in a straight-forward manner and added 
some further subtleties to incorporate higher corrections of height correlations in the fields. 
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2.1 Grassmannian Method 

In their paper, Moghimi-Araghi, Rajabpour and Rouhani |14j re-expressed the Majumdar-Dhar 
probability of clusters in the bulk, in terms of a contrived statistical system according to Berezin's 
definition of Grassmannian integrals: 

= detA c = jMjdgjexp (E^A^) = J d Q id Q j exp (E(gA^ + 0^6 j)) 
[ ' detA fdOidejexpfeOiAijOj) f dOidOj exp (E0*Ayfy) ' 

Here, A is the lattice Laplacian matrix and B is the defect matrix defined in ref. [6]. 

Evidently, since [0i0j,0kQi] = for all i, j, k and /, and hence, one can use Baker-Hausdorff- 
Campbell formula, the above expression can be interpreted as the expectation value of a field: 



Lp c = exp 



(E Wi) • ( 2 ) 



Subtleties arise in transition to the continuum limit; one can perform the transition in different 
ways which yield different results. Naively, one can find the continuum limit of TiOiBijOj in the 
exponent by expanding 9 and 9 around 9q. As an example, if we take the cluster to be a single 
height-1 site, the corresponding field obtained with this method turns out to be 

OiBijO^so ^ deBe + Beae - -ee. (3) 

As 9 and 9 are Grassmann variables, the resulted exponential can be calculated easily, since its 
Taylor series ends at quadratic terms in 9 and 9. The fields obtained in this manner are not in 
accord with results of Mahieu and Ruelle [I], though there exist some similarities. 

The second more careful method of transition to continuum is to expand the exponential in 
terms of 9 and 9 first, and expanding 9 and 9 around 9q thereafter. As an illustration, suppose 
we want to find out the field corresponding to a height-1 site in the bulk [14J. This can be sought 
through the correlation function of two height-1 fields in the bulk in Majumdar-Dhar method [6] 
and relating that to the fictitious statistical system: 

/ dStSj exp (E^iA^ + 9 i B\ i 9 j + 9 i Bj j 9 j 
P(C 1 ,C 2 ) = det(l + GB) = V " 



Jd9 i d8 j exp (E 9 l A lJ 8 j ) 

exp (OiBlfi,) exp (4) 

Here, G is inverse of the matrix A. 

If we calculate (4>Ci4>C2) using Wick theorem, we find different terms. One can classify them 
according to the number of long-range contractions between Grassmann variables of the two fields 
- this number is always even, as a contraction of a 9 of C\ to a 9 of C2 must be compensated with 
a contraction of a 8 of C2 with a 9 of C\. Having no long-rang contractions, simply reveals the 
probabilities of the clusters individually and says nothing about correlation of the two fields. The 
first relevant order is when we have two long-range contractions and contract the other Grassmann 
variables within their own clusters. So, to obtain the scaling field, on can contract all the 9 , s and 
(?'s leaving a a pair of them uncontracted. The resulting field is of the form <p = E 9iAij9j. Now, 
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expanding all 6>'s and 6>'s around 9q to the second order, one obtains the desired scaling field; for 
instance, for one site cluster with height one, we have: 

(j) So ( z ) = - 4 ^ ? 2 > ■. deBe + Beae : . (5) 

which is the same field obtained in ref. [3J. This expansion is only up to the leading term in a/r, 
with a and r being lattice spacing and and typical distance between two scaling fields. In the next 
section, we consider higher order terms and derive corrections to these scaling fields. 



3 Calculation of Higher Orders 

Obtaining the scaling fields in this way is a wearisome belabored task and is prone to be afflicted 
by human error. Thus, we developed a Mathematica [20] code by which we could handle these 
calculations fast, reliable and adjustable for different configurations. We checked the code by 
recalculating the fields assigned with simple WAC's; the final results were in accordance with those 
of Moghimi-Araghi et. al. [H] . 

Now we would like to see what fields would arise if we consider higher order terms in a/r. The 
configuration probability, P(C) and the correlation function, P(C\,C2), are the only significant 
data we have at hand to figure out the form of the corresponding scaling field. Following the 
Grassmannian method, one can consider the expression exp (9B9) as the field associated with a 
WAC. First we expand this expression in terms of 9 and 9. In Grassmannian method we contract 
all the fields except two of them, but we would like to derive several different terms depending on 
how many 9's and 9's are left uncontracted. The resulting field will have the form: 

<p c = c i + E + E ^kfiWi + ...+E AyL.^A • • • o i2n -A 2n +■■■, (6) 

Note that the above series terminates due to Grassmann nature of 9 and 9. In each term n is 
the number of uncontracted pairs, and these pairs, when contracted with similar pair of another 
field produce long-range correlations. To derive the coefficient A^ one can use a generalized 
version of graphical method used in [TJJ; yet it is possible to find these coefficients in an easier 
way. Contraction of all pairs will give us the probability of the cluster which is (4>c) = P{p) = 
det (I + BG). It is easy to see that the coefficient of Gij in P{C) is A$ and the coefficient of 

(2) 

GijGki is A\j kl . In this manner, we can find out all the needed coefficients, *4s, in (ftc easily; that is, 
firstly, we compute P(C) = det(J + G B) as a function of Gj/s. Secondly, we derive the coefficient 
of G^, for all i and j to find A^ or derive the coefficient of GijGki to find A^ kl . Then remains 
just one additional step: to expand #j's and 9^s around 9q. As an example, doing all the procedure 
explained above we arrive at the following expression for one-site height-one cluster up to 0(a 4 ): 

(j) So = a 2 4 ^ 2 7 — (9989 + B989) 

7T 

+ a 4 2 ^ 3 ~ 2 ^ (d9d 3 9 + B 3 9d9 + B9B 3 9 + B 3 9B9 - 2( ^_ ^ (b 2 + B 2 ) 9 (b 2 + B 2 ) oj 
+ a 4 32 ( 2 ~ 7r ) 89398989. (7) 

7T 
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In simplifying the original terms of the field, the evolution equations of c = —2 action was 
applied to eliminate some terms: 

dm = dBe = o. (8) 

The first line of equation (|7|) is the ordinary field derived before. The second line is the expansion of 
up to 0(a 4 ), and the third line originates from the term .A^jy- Note that there are no terms of 
the order a 3 . In fact there were some terms, but all vanish due to equation motion (equation ©). 
The fields of the order of a 4 are irrelevant under renormalization group and so could be neglected 
in the scaling limit. But since in ASM, size of lattice spacings may not be neglected in general, 
these extra terms could be relevant in some calculations. 

Though we have derived the higher corrections to the scaling field, yet it is not the end of the 
story. If we are going to take care of terms of the order of a in scaling fields, we have to treat in the 
same way with the action. The c = —2 action is obtained if we collect only the terms up to second 
order of a. So, if the size of the lattice is important, the action admits some modification. Changing 
the action will change the Green function and we already know that in calculation of higher order 
correction to correlation functions, one has to take into account the higher order terms of discrete 
Laplacian's Green function. This can be done readily in our scheme: we have the discrete version 
of the action, just we have to expand it up to order a 4 . The result is: 

A = i J (dOdO + a 2 ( K ^{d 2 6d 2 6 + B 2 9d 2 9) + ^dd0dffij\ . (9) 

Note that the new action is still quadratic and hence integrable, however it does not have conformal 
symmetry, although the off critical term is an irrelevant term and vanishes under RG. But we are 
considering the problem when the smal scale of the system can not be neglected completely, this 
forces us to keep the extra, and irrelevant, terms in the action. The last point we would like to 
mention is that the a 3 terms in the scaling field still vanish, since the application of equation of 
motion terms of order of a 5 only. 

4 Boundary Fields 

In this section we would like to derive the properties of scaling fields in ASM when a boundary is 
present at the theory. This will help us to derive finite size properties and surface critical properties 
of the model. In the context of conformal field theory, the boundary is usually taken to be the real 
line in the complex plane and the system is supposed to fill the upper half-plane. The correlation 
functions of scaling fields in this geometry reveals surface critical exponents and the finite size 
scaling properties. It is derived that under certain boundary conditions the correlation functions of 
scaling fields in this geometry is the holomorphic part of bulk correlation functions of same fields 
together with their images (which live in the lower half plane) [21]. The idea was later generalized 
to the case of logarithmic conformal field theory [221 [23] . 

Many properties of ASM in presence of such boundary has been derived. Ivashkevich calculated 
all two-point functions of all height variables along closed and open boundaries [8]. Using these 
correlation functions, Jeng identified the height variables with certain fields in c = — 2. A more 



5 



detailed treatment is given in [IS]. However, this identification is through examining the correlation 
functions. We would like to apply the Grassmaniann method to ASM with the mentioned boundary. 
Note that as this method is only for WAC's, we are not able to derive the height-2 or more one-site 
fields. 

The question we are going to address is that if we consider the model near a boundary, the 
toppling matrix, A, and the Green function G change. Since in derivation of the fields using the 
Grassmaninam method, we use Green function, the resulting field in the presence of a boundary 
maybe different from the field in the bulk. 

To begin, we would like to calculate the one point function of a WAC in the upper half boundary. 
This is given by the following expression: 



} detA / dOidOj exp 9 i A ij e j ) J dOidOj exp (£) di&ijOj) ' 

which is the same expression as in the bulk, just you should use the appropriate matrix A. Again 
B c is the matrix defined in the Majumdar-Dhar method for the specified WAC, C. So it is seen 
that the field associated with this cluster is again exp 0, H'jOj. but we should keep in mind that to 
derive the proper form of the field we use boundary Green functions which are different from the 
bulk's Green functions. Let's expand this expression and contract all the 0's and 0's to find the 
probability of the cluster. Note that as we would like to find one-point correlation function, we do 
not leave two of the variables uncontracted. The Green function of the theory with this geometry 
is obtained easily by method of images: 

G op (ri,r 2 ) = G(x,y 1 - y 2 ) - G(x,y 1 + y 2 ), 

G c \{r u f 2 ) = G(x,y x - y 2 ) + G(x,y 1 + y 2 - 1), (11) 

where G p/ci i s the Green function open/closed boundary conditions and G is the Green function 
in the bulk. Also r- L = (xi,yi) and x = x 2 — x\. 

So, each contraction of 9 variables has two terms; one, which we call it short range (SR) 
contraction, comes from the first terms in above equations and is the bulk Green function of the 
points r\ and r 2 . The other one, which we call it image- long range (ILR) contraction, is the bulk 
Green function of the point T\ with the image of the point f 2 (see figure 1). 

Take a typical expression of 9's and 9's to be fully contracted. We can classify the terms 
appearing in the result by number of ILR contractions (ILRC's), nj. If ni is set to zero, we will 
arrive at the bulk probability of the cluster C. The first correction due to presence of the boundary 
appears when nj = 1. This means that we should contract all the variables in the way we did in 
the bulk except two which we will contract later using a ILR contraction. The procedure reduces 
to the one discussed by [J3] and for example the result for the height-1 one-site cluster would be 
given by equation([S]), only you have to contract the remaining 9 and 9 using ILRC. In the scaling 
limit this gives 

<<feo(r)}o P /cl = P(l) (lT^ + -), (12) 
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Figure 1: Examples of SR-, ILR- and T-contractions near a boundary 



where y is distance of the field from the boundary and minus/plus sign corresponds to closed/open 
boundary condition which is consistent with previous results. Also if we take a typical WAC, with 
the corresponding field [U [9] 

4> s = - [A s : 8989 + 8989 : +B 1S ■ 0689 + 8989 : +iB 2S : 9090 - <9#d0 :] (13) 

it is easy to check that the second and third term of this field, do not give any contribution when 
arc ILR-contracted. Hence the one-point boundary correlation function is obtained to be 

(0s)op/ci = P(S) T ^, (14) 

just as it is indicated in [9]. 

Now we move on to calculate the field of a WAC explicitly. Suppose you have two WAC's C\ and 
C2 and would like to compute the probability of such configuration. Again we have an expression 
like (jlj). Now we can classify them depending on the number of trans-contractions (TC's), tit and 
the number of ILRC's among the inter-contractions, nj. Setting both ny and nj equal to zero, 
we will arrive at P{C\)P{C2) with P(C) being the bulk probability of the cluster C. This term 
together with the term coming from ny = and nj = 1 reveals (4>c± (n.))op/cl(0C 2 (^2))op/ch the 
disjoint boundary probability of the two clusters. To find the correlation of the two clusters, we 
should set ut nonzero. Taking tit = 2 (the smallest nonzero value for tit) and nj = the problem 
reduces to the one in the bulk: contract all 0's and 0's using bulk Green function except two, which 
are left to be contracted with the 9 variables of the other cluster's field. This means that the field 
derived in this manner is exactly the same as the one in the bulk and we do not need to care about 
the fact that the Green functions have changed. 

The first terms that contain both joint probabilities and the effect of boundary Green function, 
appears when we set tit = 2 and nj = 1. This means that you should contract all 0's 9 , s but 
four, using SRC, two of remaining will be contracted with the variables of other fields and two will 
be ILR contracted. Such field has been calculated in the previous section and we have seen it is 
of the order of a and vanishes in the scaling limit. So, if we are considering the problem in the 
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scaling limit, the Grassmaniann method says the the boundary field is just the same as bulk field. 
However, this does not mean that the effect of boundary could be neglected completely, because in 
the trans-contractions should use boundary Green functions. 
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